function [solzarray]= provideFamilySolution_matrixpoly(F,w,h,numparams,numsols)
%this function , given a camera center and a focal length and a fundamental
%matrix computes the error with respect to a fundamental matrix
if(numparams==2)
    solzarray=findFocal12(F,w,h,numparams,numsols);
elseif(numparams==3)
    solzarray=findFocalOC(F,w,h,numparams,numsols);
end
    



end


function [solzarray]=findFocal12(F,w,h,numparams,numsols)
F=F/norm(F,'fro');

focalstart=w/2;
focalend=w*2;
step=((focalend-focalstart)/(numsols-1));
solzarray=zeros(numsols,numparams);

f1=F(1,1);
f2=F(1,2);
f3=F(1,3);

f4=F(2,1);
f5=F(2,2);
f6=F(2,3);

f7=F(3,1);
f8=F(3,2);
f9=F(3,3);


if(numparams==2)
    xc=w/2;
    yc=h/2;
    
else
    display('not imeplemented yert');
end
counter=1;

xarray=focalstart:step:focalend;
for i=1:numsols
    
    foc2=xarray(i);
     p=[
    (f1^2*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)^2 - 4*f1^3*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc)))
    0
    (- 4*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)^3 - 12*f1^2*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) + 18*f1*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc)) + 2*f1*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)^2*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) + 2*f1^2*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2))
    0
    (- 12*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)^2 - 27*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))^2 + (f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)^2*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc)^2 + f1^2*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))^2 + 18*f1*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc)) + (f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(18*f5*foc2^2 + 18*f9 + 18*xc*(f7 + f1*xc + f4*yc) + 18*yc*(f8 + f2*xc + f5*yc) + 18*f3*xc + 18*f6*yc) - 12*f1*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc)^2 + 4*f1*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc))
    0
    (- 12*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))^2*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2) - 4*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc)^3 + (f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*((f1*f5*foc2^2 - f2*f4*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - f5*foc2^2*(f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f4*foc2^2*(f3 + f1*xc + f2*yc)*(f8 + f2*xc + f5*yc) + f2*foc2*(f7 + f1*xc + f4*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) - f1*foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(18*f5*foc2^2 + 18*f9 + 18*xc*(f7 + f1*xc + f4*yc) + 18*yc*(f8 + f2*xc + f5*yc) + 18*f3*xc + 18*f6*yc) + 2*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f1*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - (f3 + f1*xc + f2*yc)*(f7 + f1*xc + f4*yc) + f1*f5*foc2^2 - f2*f4*foc2^2)*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc)^2 + 2*f1*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))^2*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc))
    0
    (f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))^2*(f5*foc2^2 + f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc)^2 - 4*(f5*foc2^2*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))^3
    ];
r = roots(p');
    
solzarray(counter,1)=processRoots_focalonly(r);
solzarray(counter,2)= foc2;
    

counter=counter+1; 
end




%plot(solzarray(:,2),solzarray(:,1));


end

function [sls]=findFocalOC(F,w,h,numparams,numsols)



end

function [sl]=processRoots_focalonly(p)

s=size(p,1);
minc=1000000000;
flagset=0;
sl=0;
for i=1:s
    if(real(p(i,1))>0 && abs( imag(p(i,1)))<minc)
        sl=real(p(i,1));
        minc=abs(imag(p(i,1)));
        flagset=1;
    end
end

if(flagset==0)
    display('no solution');
end


end


% open stock.m
% clear all
% clc
% syms f1 f2 f3 f4 f5 f6 f7 f8 f9 foc1 foc2 xc yc real
% K=[foc1 0 xc ; 0 foc2 yc; 0 0 1]
% F=[f1 f2 f3 ; f4 f5 f6 ; f7 f8 f9 ]
% G=(K')*F*K
% v=poly(G'*G)
% a=1
% b=(- f1*foc1^2 - f5*foc2^2 - f9 - xc*(f7 + f1*xc + f4*yc) - yc*(f8 + f2*xc + f5*yc) - f3*xc - f6*yc)
% c=((f1*foc1^2 + f5*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) - foc1*(f7 + f1*xc + f4*yc)*(f3*foc1 + f1*foc1*xc + f2*foc1*yc) - foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) + f1*f5*foc1^2*foc2^2 - f2*f4*foc1^2*foc2^2)
% d= - (f1*f5*foc1^2*foc2^2 - f2*f4*foc1^2*foc2^2)*(f9 + xc*(f7 + f1*xc + f4*yc) + yc*(f8 + f2*xc + f5*yc) + f3*xc + f6*yc) + (foc1*(f7 + f1*xc + f4*yc)*(f3*foc1 + f1*foc1*xc + f2*foc1*yc) + foc2*(f8 + f2*xc + f5*yc)*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f1*foc1^2 + f5*foc2^2) - foc1*(f1*foc1^2*(f3*foc1 + f1*foc1*xc + f2*foc1*yc) + f2*foc1*foc2*(f6*foc2 + f4*foc2*xc + f5*foc2*yc))*(f7 + f1*xc + f4*yc) - foc2*(f5*foc2^2*(f6*foc2 + f4*foc2*xc + f5*foc2*yc) + f4*foc1*foc2*(f3*foc1 + f1*foc1*xc + f2*foc1*yc))*(f8 + f2*xc + f5*yc)
% delta=18*a*b*c*d-(4*b^3*d)+(b^2*c^2)-(4*a*c^3)-(27*a^2*d^2)
% dsim=simplify(delta)
% RS = collect(dsim,foc1)